/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 1999, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** ****************************************************************** */
                                                                       
// $Revision: 1.6 $
// $Date: 2009-11-02 22:23:58 $
// $Source: /usr/local/cvs/OpenSees/EXAMPLES/Example1/main.cpp,v $


// File: ~/model/main.C
//
// Written: fmk 08/99
//
// Purpose: this file contains a C++ main procedure to perform the analysis
// of example1 (found in most documents). In the main() procedure:
//      1) each object of the domain, i.e. Nodes, Elements, Constraints,
//         and LoadPattern objects are created and then added to the Domain.
//      2) the components of the analysis object are constructed and then
//         the Analysis object is created.
//      3) the analysis is performed.
//      4) the results are printed - here the contents of Domain and end of
//         the analysis operation.

// standard C++ includes

#include <stdlib.h>

#include <OPS_Globals.h>
#include <StandardStream.h>

#include <ArrayOfTaggedObjects.h>

// includes for the domain classes
#include <Domain.h>
#include <Node.h>
#include <Truss.h>
#include <ElasticTimoshenkoBeam2d.h>
#include <ElasticBeam2d.h>
#include <ElasticMaterial.h>
#include <SP_Constraint.h>
#include <LoadPattern.h>
#include <LinearSeries.h>
#include <NodalLoad.h>
#include <CrdTransf.h>
#include <PDeltaCrdTransf2d.h>

// includes for the analysis classes
#include <StaticAnalysis.h>
#include <AnalysisModel.h>
//#include <Linear.h>
#include <KrylovNewton.h>
#include <LineSearch.h>
#include <NewtonLineSearch.h>
#include <InitialInterpolatedLineSearch.h>
#include <TransformationConstraintHandler.h>
#include <DOF_Numberer.h>
#include <RCM.h>
#include <LoadControl.h>
#include <BandSPDLinSOE.h>
#include <BandSPDLinLapackSolver.h>
#include <SymSparseLinSOE.h>
#include <SymSparseLinSolver.h>
#include <CTestEnergyIncr.h>
#include <CTestNormDispIncr.h>
#include <ID.h>
#include <OPS_Stream.h>
#include <NodeRecorder.h>
#include <Recorder.h>

// init the global variabled defined in OPS_Globals.h
StandardStream sserr;
OPS_Stream *opserrPtr = &sserr;




// main routine
int main(int argc, char **argv)
{
    //
    //  now create a domain and a modelbuilder
    //  and build the model
    //

    Domain *theDomain = new Domain();
   
    // create the nodes using constructor:
    //          Node(tag, ndof, crd1, crd2)
    // and then add them to the domain
    Vector Trans1(2);
    Trans1(0) = 0;
    Trans1(1) = 0;
    Vector Trans2(2);
    Trans2(0) = 0;
    Trans2(1) = 0;
    CrdTransf *PDelta1 = new PDeltaCrdTransf2d(1);

    Node *node1 = new Node(1, 3,   0.0,  0.0);
    Node *node2 = new Node(2, 3, 1440.0,  0.0);
    Node *node3 = new Node(3, 3, 1680.0,  0.0);
    Node *node4 = new Node(4, 3,  720.0, 960.0);
    //Node *node5 = new Node(5, 3, 0.0, -100.0);
    theDomain->addNode(node1);
    theDomain->addNode(node2);
    theDomain->addNode(node3);
    theDomain->addNode(node4);
    //theDomain->addNode(node5);
   
    // create an elastic material using constriuctor:  
    //          ElasticMaterialModel(tag, E)

    UniaxialMaterial *theMaterial = new ElasticMaterial(1, 3000);
   
    // create the truss elements using constructor:
    //          Truss(tag, dim, nd1, nd2, Material &,A)
    // and then add them to the domain
   
    Truss *Beam1 = new Truss(1, 2, 1, 4, *theMaterial, 11600);
    //Truss2 *truss2 = new Truss2(2, 2, 2, 4, *theMaterial,  5.0);
    //Truss2 *truss3 = new Truss2(3, 2, 3, 4, *theMaterial,  5.0);
    //Truss *Beam1 = new Truss(4, 2, 5, 1, *theMaterial,  5.0);

    //ElasticTimoshenkoBeam2d *Beam1 = new ElasticTimoshenkoBeam2d(4, 5, 1, 3000, 2000, 8, 0.08, 7, *PDelta1);
    //ElasticTimoshenkoBeam2d *Beam1 = new ElasticTimoshenkoBeam2d(1, 1, 4, 206000, 206000, 11600, 162786666, 10000, *PDelta1);
    //ElasticTimoshenkoBeam2d *Beam2 = new ElasticTimoshenkoBeam2d(2, 2, 4, 206000, 206000, 11600, 162786666, 10000, *PDelta1);
    //ElasticTimoshenkoBeam2d *Beam3 = new ElasticTimoshenkoBeam2d(3, 3, 4, 206000, 206000, 11600, 162786666, 10000, *PDelta1);

    //ElasticBeam2d *Beam1 = new ElasticBeam2d(1, 11600, 206000, 162786666, 1, 4, *PDelta1);
    ElasticBeam2d *Beam2 = new ElasticBeam2d(2, 11600, 206000, 162786666, 2, 4, *PDelta1);
    ElasticBeam2d *Beam3 = new ElasticBeam2d(3, 11600, 206000, 162786666, 3, 4, *PDelta1);

    theDomain->addElement(Beam1);
    theDomain->addElement(Beam2);
    theDomain->addElement(Beam3);
    //theDomain->addElement(Beam1);
   
    // create the single-point constraint objects using constructor:
    //          SP_Constraint(tag, nodeTag, dofID, value)
    // and then add them to the domain
    SP_Constraint *sp4 = new SP_Constraint(2, 0, 1);
    SP_Constraint *sp5 = new SP_Constraint(2, 1, 1);
    //SP_Constraint *sp6 = new SP_Constraint(2, 2, 1);
    SP_Constraint *sp7 = new SP_Constraint(3, 0, 1);
    SP_Constraint *sp8 = new SP_Constraint(3, 1, 1);
    //SP_Constraint *sp9 = new SP_Constraint(3, 2, 1);
    SP_Constraint *sp1 = new SP_Constraint(1, 0, 1);
    SP_Constraint *sp2 = new SP_Constraint(1, 1, 1);
    SP_Constraint *sp3 = new SP_Constraint(1, 2, 1);
    theDomain->addSP_Constraint(sp1);
    theDomain->addSP_Constraint(sp2);
    theDomain->addSP_Constraint(sp3);
    theDomain->addSP_Constraint(sp4);    
    theDomain->addSP_Constraint(sp5);    
    //theDomain->addSP_Constraint(sp6);
    theDomain->addSP_Constraint(sp7);
    theDomain->addSP_Constraint(sp8);
    //theDomain->addSP_Constraint(sp9);

    // construct a linear time series object using constructor:
    //          LinearSeries()
   
    TimeSeries *theSeries = new LinearSeries();
   
    // construct a load pattren using constructor:
    //          LoadPattern(tag)
    // and then set it's TimeSeries and add it to the domain
   
    LoadPattern *theLoadPattern = new LoadPattern(1);
    theLoadPattern->setTimeSeries(theSeries);
    theDomain->addLoadPattern(theLoadPattern);
   
    // construct a nodal load using constructor:
    //          NodalLoad(tag, nodeID, Vector &)
    // first construct a Vector of size 2 and set the values NOTE C INDEXING
    // then construct the load and add it to the domain
   
    
    Vector theLoadValues(3);
    theLoadValues(0) = 10000.0;
    theLoadValues(1) = -5000.0;
    theLoadValues(2) = 0.0;
    NodalLoad *theLoad = new NodalLoad(1, 4, theLoadValues);
    theDomain->addNodalLoad(theLoad, 1);

    ID vdofs(6);
    ID vnodes(4);
    vnodes(0) = 1;
    vnodes(1) = 2;
    vnodes(2) = 3;
    vnodes(3) = 4;
    vdofs(0) = 0;
    vdofs(1) = 1;
    //OPS_Stream *Unknown = new OPS_Stream(1,2);
    //Recorder *DispRecord = new NodeRecorder(vdofs, vnodes, 0, "disp", *theDomain, *Unknown);

    // create an Analysis object to perform a static analysis of the model
    //  - constructs:
    //    AnalysisModel of type AnalysisModel,
    //    EquiSolnAlgo of type Linear
    //    StaticIntegrator of type LoadControl
    //    ConstraintHandler of type Penalty
    //    DOF_Numberer which uses RCM
    //    LinearSOE of type Band SPD
    // and then the StaticAnalysis object
   
    AnalysisModel     *theModel = new AnalysisModel();
    StaticIntegrator  *theIntegrator = new LoadControl(0.1, 10, 0.1, 0.1);
    ConstraintHandler *theHandler = new TransformationConstraintHandler();
    RCM               *theRCM = new RCM();

    DOF_Numberer      *theNumberer = new DOF_Numberer(*theRCM);    
    //BandSPDLinSolver  *theSolver = new BandSPDLinLapackSolver();
    //LinearSOE         *theSOE = new BandSPDLinSOE(*theSolver);
    SymSparseLinSolver  *theSolver = new SymSparseLinSolver();
    SymSparseLinSOE *theSOE = new SymSparseLinSOE(*theSolver, 1);
    ConvergenceTest   *theTest = new CTestEnergyIncr(0.000001, 200, 1, 2, 1);

    LineSearch *theSearch = new InitialInterpolatedLineSearch();
    EquiSolnAlgo      *theSolnAlgo = new NewtonLineSearch(*theTest, theSearch);
    StaticAnalysis    theAnalysis(*theDomain, *theHandler, *theNumberer, *theModel, *theSolnAlgo, *theSOE, *theIntegrator, theTest);

    // perform the analysis & print out the results for the domain
    int numSteps = 50;
    int nstep = 0;
    for (nstep = 0; nstep < numSteps; nstep++)
    {
        theAnalysis.analyze(1);
    }
    ConvergenceTest *theTest2 = new CTestNormDispIncr(0.000001, 200, 1, 2, 1);
    StaticAnalysis theAnalysis2(*theDomain, *theHandler, *theNumberer, *theModel, *theSolnAlgo, *theSOE, *theIntegrator, theTest2);

    for (nstep = 0; nstep < numSteps; nstep++)
    {
        theAnalysis2.analyze(1);
    }
    opserr << *theDomain;
    exit(0);
}
